Phonon-induced spin-spin interactions in diamond nanostructures: 
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We propose and analyze a novel mechanism for long-range spin-spin interactions in diamond 
nanostructures. The interactions between electronic spins, associated with nitrogen-vacancy centers 
in diamond, are mediated by their coupling via strain to the vibrational mode of a diamond me- 
chanical nanoresonator. This coupling results in phonon-mediated effective spin-spin interactions 
that can be used to generate squeezed states of a spin ensemble. We show that spin dephasing 
and relaxation can be largely suppressed, allowing for substantial spin squeezing under realistic 
experimental conditions. Our approach has implications for spin-ensemble magnetometry, as well 
as phonon-mediated quantum information processing with spin qubits. 
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Electronic spins associated with nitrogen-vacancy 
(NV) centers in diamond exhibit long coherence times 
and optical addressability, motivating extensive research 
on NV-based quantum information and sensing applica- 
tions. Recent experiments have demonstrated coupling 
of NV electronic spins to nuclear spins [TJ [2], entangle- 
ment with photons [3] , as well as single spin |U [5] and 
ensemble [6l [7] magnetometry. An outstanding challenge 
is the realization of controlled interactions between sev- 
eral NV centers, required for quantum gates or to gener- 
ate entangled spin states for quantum-enhanced sensing. 
One approach toward this goal is to couple NV centers 
to a resonant optical [SJ [S] or mechanical [Tt)HT2"] mode; 
this is particularly appealing in light of rapid progress in 
the fabrication of diamond nanostructures with improved 
optical and mechanical properties |13H17j . 

In this Letter, we describe a new approach for effec- 
tive spin-spin interactions between NV centers based on 
strain-induced coupling to a vibrational mode of a dia- 
mond resonator. We consider an ensemble of NV cen- 
ters embedded in a single crystal diamond nanobeam, as 
depicted in Fig. [IJl When the beam flexes, it strains 
the diamond lattice which in turn couples directly to 
the spin triplet states in the NV electronic ground state 
[TBI ITU] . For a thin beam of length L ~ 1 /im, this 
strain-induced spin-phonon coupling can allow for co- 
herent effective spin-spin interactions mediated by vir- 
tual phonons. Based on these effective interactions, we 
explore the possibility to generate spin squeezing of an 
NV ensemble embedded in the nanobeam. We account 
for spin dephasing and mechanical dissipation, and de- 
scribe how spin echo techniques and mechanical driv- 
ing can be used to suppress the dominant decoherence 
processes while preserving the coherent spin-spin inter- 
actions. Using these techniques we find that significant 



spin squeezing can be achieved with realistic experimen- 
tal parameters. Our results have implications for NV en- 
semble magnetometry, and provide a new route toward 
controlled long-range spin-spin interactions. 

Model. — The electronic ground state of the negatively 
charged NV center is a spin S — 1 triplet with spin states 
labeled by \m s = 0, ±1) as shown in Fig. [T|d. In the pres- 
ence of external electric and magnetic fields E and B, 
the Hamiltonian for a single NV is (H = 1) [T5] 

#nv =(Do+d\\E z )S 2 z +^ B g s S-B 

— d± [E x (S x S y + S y S x ) + E y {S x ~ Sy)] , 

where Do/2ir ~ 2.88 GHz is the zero field splitting, 
g s ~ 2, fis is the Bohr magneton, and rfn (d±) is the 
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FIG. 1. (a) All-diamond doubly clamped mechanical res- 
onator with an ensemble of embedded NV centers, (b) Spin 
triplet states of the NV electronic ground state. Local perpen- 
dicular strain induced by beam bending mixes the |±1) states, 
(c) A collection spins in the two- level subspace {| + 1} , | — 1)} is 
off-resonantly coupled to a common mechanical mode giving 
rise to effective spin-spin interactions, (d) Squeezing of the 
spin uncertainty distribution of an NV ensemble. 
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ground state electric dipole moment in the direction par- 
allel (perpendicular) to the NV axis [20l [21] . 

Motion of the diamond nanoresonator changes the lo- 
cal strain at the position of the NV center, which results 
in an effective, strain-induced electric field [TS]. We are 
interested in the near-resonant coupling of a single reso- 
nant mode of the nanobeam to the |±1) transition of the 
NV, with Zeeman splitting = g s /J,BB z /h, as shown 
in Fig. [T]d,c. The perpendicular component of strain 
E± mixes the |±1) states. For small beam displace- 
ments, the strain is linear in its position and we write 
E± = Eo(a + cl 1 ), where a is the destruction operator of 
the resonant mechanical mode of frequency u> m , and Eq 
is the perpendicular strain resulting from the zero point 
motion of the beam. We note that the parallel component 
of strain shifts both states |±1) relative to |0) [22]; how- 
ever, with near-resonant coupling A = A^ — w m -C Do 
and preparation in the |±1) subspace, the state |0) re- 
mains unpopulated and parallel strain plays no role in 
what follows. Within this two-level subspace, the in- 
teraction of each NV is Hi = g (afa + cra^), where 
u i = |±1) ■ (=F 1 1 is the Pauli operator of the zth NV cen- 
ter and g is the single phonon coupling strength. For 
many NV centers we introduce collective spin operators, 



J* = lE,|i)i(i| 



T| and J± = J x ± iJ y 



Ei °i i wmcn satisfy the usual angular momentum com- 
mutation relations. The total system Hamiltonian can 
then be written as 



H = uj m <va + AbJ z + g (arJ— + aJ+) 



(2) 



which describes a Tavis-Cummings type interaction be- 
tween an ensemble of spins and a single mechanical mode. 
In Eq. (|]) we have assumed uniform coupling of each spin 
to the mechanical mode for simplicity. In general the 
coupling may be nonuniform and we discuss this further 
below. 

To estimate the coupling strength g, we calculate the 
strain for a given mechanical mode and use the experi- 
mentally obtained stress coupling of 0.03 Hz Pa" 1 in the 
NV ground state [23] [24] . We take a doubly clamped 
diamond beam (see Fig. [T^,) with dimensions L 3> w,h 
such that Euler-Bernoulli thin beam elasticity theory is 
valid [S3] . For NV centers located near the surface of the 
beam we obtain 1231 
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GHz, 



(3) 



where p is the mass density and E is the Young's 
modulus of diamond. For a beam of dimensions 
(L,w,h) = (1, 0.1, 0.1) /im we obtain a vibrational fre- 
quency uj m /2iT ~ 1 GHz and coupling g/2n ~ 1 kHz. 
While this is smaller than the strain coupling g e /2ir s» 10 
MHz expected for electronic excited states of defect cen- 
ters HZ] or quantum dots [3H] , we benefit from the 
much longer spin coherence time Ti in the ground state. 



An important figure of merit is the single spin cooper- 
ativity i] — ' ; where 7 = us m /Q is the mechanical 
damping rate and nth = (e h ^™l kBT — l)" 1 is the equilib- 
rium phonon occupation number at temperature T. As- 
suming Q = 10 6 , T2 = 10 ms and T = 4 K, we obtain a 
single spin cooperativity of n ~ 0.8. This can be further 
increased by reducing the dimensions of the nanobeam 
and operating at lower temperatures. 

Spin squeezing. — In the dispersive regime, g -C A = 
As — u> m , virtual excitations of the mechanical mode 
result in effective interactions between the otherwise de- 
coupled spins. In this limit, H can be approximately 
diagonalized by the transformation e R He~ R with R — 
^ (a^J_ — aJ-f). To order (g/A) 2 this yields an effec- 
tive Hamiltonian, 

H cff = iu m a)a + (A B + Aa f a) J z + ^J+J-, (4) 

where A = 2g 2 / A is the phonon-mediated spin-spin cou- 
pling strength. Rewriting J + J_ = J 2 — J 2 + J z , and 
provided the total angular momentum J is conserved, 
we obtain a term cx J 2 corresponding to the one-axis 
twisting Hamiltonian [29] . 

To generate a spin squeezed state, we initialize the en- 
semble in a coherent spin state (CSS) \i/}q) along the 
x axis of the collective Bloch sphere. The CSS satis- 
fies J x = J \ipo) and has equal transverse variances, 
(Jy) = \J£) — J/2- This can be achieved using opti- 
cal pumping and global rotations of the spins with mi- 
crowave fields [30]. The squeezing term cx J 2 describes 
a precession of the collective spin about the z axis at a 
rate proportional to J z , resulting in a shearing of the un- 
certainty distribution and a reduced spin variance in one 
direction as shown in Fig. [l]i. This is quantified by the 
squeezing parameter (3TJ [32] , 
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where (AJ^ in ) = \ [V+ - ^JV 2 + V 2 z j is the mini- 
mum spin uncertainty with V± — (J 2 ± J 2 ) and V yz — 
(JyJz + JzJy) /2. The preparation of a spin squeezed 
state, characterized by £ 2 < 1, has direct implications for 
NV ensemble magnetometry applications, since it would 
enable magnetic field sensing with a precision below the 
projection noise limit [31 j - 

We now consider spin squeezing in the presence of real- 
istic decoherence. In addition to the coherent dynamics 
described by H c g , we account for mechanical dissipation 
and spin dephasing using a master equation [23] 



P 



^J 2 + (AB + Xa^a)J z ,p\+^Y, V ^P 

(6) 
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where T>[c]p — cpc^ — \ (c^cp + pc'c) and the single spin 
dephasing T^ 1 is assumed to be Markovian for simplic- 
ity (see below). Note that we absorbed a shift of A/2 
into Ab, and ignored single spin relaxation as T\ can 
be several minutes at low temperatures [33 . The sec- 
ond line describes collective spin relaxation induced by 
mechanical dissipation, with T 7 = -fg 2 /A 2 . Finally, the 
phonon number n = a^a shifts the spin frequency, acting 
as an effective fluctuating magnetic field which leads to 
additional dephasing. 

Let us for the moment ignore fluctuations of the 
phonon number n\ we address these in detail below. 
Starting from the CSS we plot the squeezing pa- 

rameter in Fig. [2^ for an ensemble of N = 100 spins 
and several values of n t hj in the presence of dephasing 
2TT 1 and collective relaxation T 7 . Here we calculated £ 2 
by solving Eq. Q using an approximate numerical ap- 
proach treating T 7 and T 2 separately, and verified that 
the approximation agrees with exact results for small N 
[23j . To estimate the minimum squeezing, we linearize 
the equations of motion for the averages and variances of 
the collective spin operators (see dashed lines in Fig. [2^,) . 
From these linearized equations, in the limits of interest, 
J ^> 1, fj t h 3> 1 and to leading order in both sources of 
decoherence, we obtain approximately 



2 _ 4r 7 n th 
4 ~ JXH T 2 ' 



(7) 



Optimizing t and the detuning A, we obtain the optimal 
squeezing parameter, 
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opt 



(8) 



at time t opt = T^/yJjj, similar to results for atomic sys- 
tems [34 36 . Note that for non-Markovian dephasing, 
the scaling is even more favorable [37]. In Fig. we 
plot the scaling of the squeezing parameter with J for 
small but finite decoherence, and find agreement with 
Eq. For comparison we also plot the unitary result 
in the absence of decoherence, scaling as £ 2 pt ~ J~ 2 / 3 
and limited by the Bloch sphere curvature [2"§] . 

Phonon number fluctuations. — In Eq. Q we see that 
the phonon number n — a^a couples to J z , leading to ad- 
ditional dephasing due to thermal number fluctuations. 
On the other hand, this same coupling can also lead to 
additional spin squeezing from cavity feedback, by driv- 
ing the mechanical mode [34Tl36j . In the following, we 
consider a twofold approach to mitigate thermal spin de- 
phasing while preserving the optimal squeezing. First, 
we apply a sequence of global spin echo control pulses 
to suppress dephasing from low-frequency thermal fluc- 
tuations. This also extends the effective coherence time 
T 2 of single NV spins j3U]. Second, we consider driving 
the mechanical mode, and identify conditions when this 
results in a net improvement of the squeezing. 
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FIG. 2. (a) Spin squeezing parameter versus scaled precession 
time with N = 100 spins. Solid blue lines show the calculated 
squeezing parameter for T2 = 10 ms and values of n t h as 
shown. For each curve, we optimized the detuning A to obtain 
the optimal squeezing. Blue dashed lines are calculated from 
the linearized equations for the spin operator averages. Black 
solid (dashed) line shows exact (linearized) unitary squeezing, 
(b) Optimal squeezing versus number of spins. Lower (upper) 
red line shows power law fit for n t h = 1 (10) and 22 = 1 
(0.01) s. The detuning A is optimized for each point. Other 
parameters in both plots are uj m /2n = 1 GHz, g/2n = 1 kHz, 
Q = 10 6 . 



To simultaneously account for thermal dephasing, 
driven feedback squeezing, and spin control pulse se- 
quences, we write the interaction term in Eq. Q in the 
so-called "toggling frame" [3"5] , 



fltat (t) = XJ z f(t)Sn(t). 



(9) 



The function f(t) periodically inverts the sign of the in- 
teraction as shown in the inset of Fig. [3^,, describing the 
inversion of the collective spin J z — ¥ —J z with each it 
pulse of the spin echo sequence. Phonon number fluctu- 
ations are described by 8n(t) = n(t) — n, where n is the 
mean phonon number and we have omitted an average 
frequency shift proportional to n in Eq. ([9]). The num- 
ber fluctuation spectrum S n (ui) — f dte (Sn(t)Sn(0)) 
is plotted in Fig. [3^ for a driven oscillator coupled to a 
thermal bath [2"3"j . 

We calculate the required spin moments within the 
Gaussian approximation for phonon number fluctuations, 
and obtain [231 



(J+W)=e- x (e- v(J =- 1/2) J + (0)), 



(10) 



and similar results for (J+(t)\ and (J + (t)J z (t)) . In 
Eq. ( 10 ) the dephasing parameter \ an d effective squeez- 



ing via p are given by 



2 , (1ujF(ut) - 
X = X I — 2 S n {w), 



p = X 



2tt uj 2 



2 / duj K(ujt) 
2n oj 2 



(11) 
(12) 



where S„(lj) — (S n (u>) + S n (— uj)) /2 and A n (u>) = 
(S n (w) - S n (-u)) /2. The filter function F{lut) = 
^- If dte lu}t f {t)\ describes the effect of the spin echo 
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pulse sequence with time r between 7r pulses 
The function K(wt) plays the analogous role for the ef- 
fective squeezing described by fi, and is related to F by 
a Kramers-Kronig relation [23] • We plot K and F for a 
sequence of M — 4 pulses in Fig. [3^,. 

Discussion. — We now consider the impact of thermal 
fluctuations on the achievable squeezing. The noise spec- 
trum S n (uj) — 27fith(^th + l)/(^ 2 + 7 2 ) is symmetric 
around ui = 0. Without spin echo control pulses, this 
low frequency noise results in nonexponcntial decay of 
the spin coherence, xo(t) = § A 2 f?.^ h t 2 (with n t h ~> 1), 
familiar from single qubit decoherence [30] S2] . The in- 
homogeneous thermal dephasing time is T 2 * ~ v^/Anth, 
severely limiting the possibility of spin squeezing. In par- 
ticular, at time t = t opt we find that squeezing is prohib- 
ited when n t h > \/J |23j . However, one can overcome 
this low frequency thermal noise using spin echo. By ap- 
plying a sequence of M equally spaced global 7r-pulses 
to the spins during precession of total time t, we obtain 
Xth ~ A 2 7n 2 h i 3 /M 2 , suggesting that thermal dephasing 
can be made negligible relative to both T 7 and 2TT 1 . For 
a sufficiently large number of pulses, M ^> fithVl^, we 
recover the optimal squeezing in Eqs. ^ and ([8]). 

Adding a mechanical drive can further enhance squeez- 
ing via feedback; however, it also increases phonon num- 
ber fluctuations, contributing to additional dephasing. 
We consider a detuned external drive of frequency cjdr = 
^Vn + leading to two additional peaks in S n (tu) at 
uj = ±5, as shown in Fig. [3ji. The area under the left 
[right] peak scales as n^nth [^dr(^th + 1)]> where nd r 
is the mean phonon number due to the drive at zero 
temperature. The symmetric and antisymmetric parts 
of this noise contribute to dephasing and squeezing as 



described by Eqs. (Ill and (12 1. Choosing the inter- 
val t/M — 2tt/S between tt pulses, we obtain additional 
dephasing Xdv — (f ) ndrnthjt and effective squeezing 
with /i ~ -g-rifct. In the limit n^ x 3> n t h, the effects of 
the drive dominate over Xth and T 7 and we recover the 
ideal scaling given in Eq. j8j, even with a small number 
of echo pulses. This is shown in Fig. [3}j,c where we see 
that the optimal squeezing improves with increasing n^- 
for a fixed number of pulses M = 4. 

Finally, we discuss our assumption of uniform coupling 
strength g in Eq. |2]). This is an important practical is- 
sue, as we expect the coupling to individual spins to be 
inhomogeneous in experiment due to the spatial variation 
of strain in the beam. Nonetheless, even with nonuni- 
form coupling, we still obtain squeezing of a collective 
spin with a reduced effective total spin J c g < J, pro- 
vided J > 1. First, we note that inhomogeneous mag- 
netic fields resulting in nonuniform detuning are com- 
pensated by spin echo. Second, for a distribution of cou- 
pling strengths gi, the effective length of the collective 
spin is JA Si/vE; .9 2 f° r the direct squeezing term, and 
Si 9i I V^2i 9i f° r feedback squeezing with a mechanical 
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FIG. 3. (a) Number fluctuation spectrum of thermal driven 
oscillator. Center (blue) peak is purely thermal while side 
(green) peaks are due to detuned drive. Solid (dashed) purple 
line shows filter function F (K) for M — 4 pulses. Inset: 
corresponding function f(t) for M — 4. (b) Solid green curves 
show squeezing parameter versus precession time for n t h = 10 
and n<ir = 10 3 ,5 x 10 4 , 10 6 (top to bottom). Dashed black 
line shows unitary squeezing, (c) Minimum squeezing versus 
drive strength for nth = 50, 10 (top to bottom). Symbols 
mark corresponding points with (b). Dashed black line shows 
unitary squeezing. Parameters in (b) and (c) are M = 4, 
g/2-K = 1 kHz, T 2 = 10 ms, TV = 100, cj m /27r = 1 GHz, 
Q = 10 6 . 



drive. Similar conclusions were reached in atomic and 
nuclear systems [54Tf3o1 14"3] . In the case of direct squeez- 
ing, it is important that the sign of the g^s is the same 
to avoid cancellation; this is automatically achieved by 
using NV centers implanted on the top of the beam. For 
beam dimensions (1, 0.1, 0.1) /mi analyzed above, we esti- 
mate that N ~ 200 NV centers can be embedded without 
being perturbed by direct magnetic dipole-dipole interac- 
tions. A reduction of the effective spin length by factor 
~ 2 still leaves N c g ~ 100, sufficient to observe spin 
squeezing. 

Conclusions. — We have shown that direct spin-phonon 
coupling in diamond can be used to prepare spin squeezed 
states of an NV ensemble embedded in a nanoresonator, 
even in the presence of dephasing and mechanical dissi- 
pation. With further reductions in temperature, beam 
dimensions, and spin decoherence rates, the regime of 
large single spin cooperativity 77 3> 1 could be achieved. 
This would allow for coherent phonon-mediated interac- 
tions and quantum gates between two spins embedded 



in the same resonator via H\ r 



X (a^ (T~ 2 



h.c. ) , and 



coupling over larger distances by phononic channels [26] . 
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SUPPLEMENTAL INFORMATION 
Coupling strength 

We assume that the NV axis is aligned with both the 
magnetic field and with the direction of beam deflection, 
so that the longitudinal strain due to deflection is entirely 
perpendicular to the NV axis. From experiment [T], the 
splitting of the |±1) states with stress is ~ 0.03 Hz/Pa. 
We convert this into the deformation potential coupling 
frequency, T E,/2nh — 36 GHz, using the Young's modulus 
of diamond, E = 1200 GPa. Next we calculate the strain 
at the NV center using elasticity theory. The equation 
for the bending mode of a thin beam is 



(z,t) = -EId*<f>{z,t) 



(13) 



where <j> is the transverse displacement in the x direction 
and z is along the beam. Here p is the mass density, A is 
the cross sectional area, and I = wt 3 /12 is the moment of 
inertia. The solutions are of the form <fi(z,t) — u(z)e~ lu}t 
where 



u{z) 



1 



cos kz — cosh kz 

(cos kL — cosh kL) 
(sin kL — sinh kL) 



(sin kz — sinh kz) 



(14) 



which satisfies the boundary conditions u(0) = u'(0) = 
u(L) = u'(L) = for a doubly clamped beam. The 
allowed wavenumbers k n are given by the solutions of 
cos kL cosh kL — 1, and the corresponding eigenfrequen- 
cies are 



uj n = k 2 n 




(15) 



We normalize the modefunction of the fundamental mode 
uq(z) by setting the free energy stored in the beam to the 
zero point energy, 



w = Iei 



dz 



cPuo 
dz 2 



Hujo 
~2~' 



(16) 



Integrating by parts we obtain the normalization condi- 
tion, 



dzuQ 



pAu 



(17) 



If the NV center lies at the midpoint along the beam, 
z = L/2, and at a distance r from the neutral axis of 
the beam, the strain due to the zero point motion of the 
fundamental mode is 



(q = -r d^u (L/2) ~ r 4 



27 



.2r 



pAuj L 5 / 2 t 



L 3 w\/Ep' 
(18) 
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where we Eq. (15) and fco ~ 4.73/L for the fundamental 



mode. The coupling strength is given by the deformation 
potential and the strain due to zero point motion, 



2tt 2wh 



e ~ 180 GHz 



2r 



L 3 W«/Ep" 



(19) 



For an NV near the surface of the beam, r ~ t/2 and 
we obtain Eq. (4) of the main text. 



Effective squeezing Hamiltonian from spin-phonon 
coupling 

In this section we provide additional details on deriving 
H e ff in Eq. (4) from the original i?NV in Eq. (1). Assum- 
ing that the magnetic field is aligned along the NV axis, 



B = Bz, and defining E± = 
we can rewrite -£?nv as (h - 



E x HE y and S± 
= 1) 



S x ziziS v 



#nv = (A) + d\\E s )Sl + g s p B BS z - 



d± 



[E + Sl 



E_S 2 _ 
(20) 



where E is the effective electric field due to strain. We 
quantize the perpendicular strain field, E + = E$a and 
= E$a\ where Eq is the strain due to the zero point 
motion of the resonant mode. Next we focus on the two- 
level subspace {|1) , | — 1)} only, and assume that transi- 
tions to state |0) are not allowed due to the large zero 
field splitting Do- For the ith spin we write Pauli opera- 
tors of = |±1>, (Tl| and af = (1| - |-1), ; (-1|, and 
within this two-level subspace the interaction for a single 
NV is 

Hi = ^pof + 9 {<rfa + oVr) + w m ata, (21) 

where g = —d±E , A B = 2g s p B B is the energy between 
|±1), and we included the mechanical oscillator of fre- 
quency uj m . Summing Eq. (21 1 for many NVs coupled to 



the same mode with uniform coupling strength we obtain 

H = A B J Z + g (at J- + aJ + ) + w m a f a, (22) 

which is Eq. (2) of the main text. To obtain Eq. (4), we 
first rewrite H in the rotating frame at the mechanical 
frequency cj m , 



H = AJ Z + g (a f J_ + aJ + 



(23) 



where A = As — u) m . Next we apply the transforma- 
tion e R He~ R , with R = ^ (a' J- — a>J+), and to order 
(g/A) 2 we obtain 



eflf 



AJ Z + ^ (J+J- + 2ataJ z ) 



(24) 



Transforming back to the nonrotating frame yields 
Eq. (4) of the main text. 



Individual spin dephasing and phonon-induced 
relaxation 



Individual spin dephasing from intrinsic T2 

Each NV spin experiences intrinsic decoherence in the 
absence of the mechanical mode. Individual relaxation 
(Ti) processes are due to lattice phonons; at low temper- 
ature T\ can be ~ 100 s and we ignore it [2]. However, 
we include intrinsic single spin dephasing, which arises 
from magnetic noise of 13 C nuclear spins in the diamond 
lattice. In practice, single spin dephasing may be non- 
exponential [21 3] , but for simplicity we approximate the 
effect of single spin dephasing by an effective Markovian 
master equation with dephasing rate 



1 

2T 2 



(25) 



Collective phonon-induced spin relaxation 

The transformation R used to obtain the effective 
squeezing Hamiltonian H c g also introduces a relaxation 
channel for the collective spin by admixing phonon and 
spin degrees of freedom. Mechanical dissipation is de- 
scribed by the master equation for the system density 
matrix p, 



f> =7(^th + 1) \apat — \ (a'ap + pa' a)] 



-frith [at pa 



aa' p - 



paat )] 



(26) 



Transforming a and using the transformation R in the 
main text, we obtain effective spin relaxation terms in 
the master equation, 



7 (n th + 1) [J-pJ+ - \ {J+J-P + pJ+J-)] 
T 7 n th [J+pJ- - I (J-J + p + pJ-J+)] 



(27) 



where 1% = (£) 7. From Eqs. (|25| and \27\ we calculate 



the equations for spin averages and variances accounting 
for individual dephasing and collective relaxation using 
dt (A) = tr{Ap}. 



Squeezing estimate from linearized equations for spin 
averages and variances 

Here we sketch the derivation of the estimated optimal 
squeezing given in Eq. (8) in the main text. In order 
to treat the squeezing Hamiltonian, collective relaxation 
and spin dephasing on equal footing, we linearize the 
equations for the spin averages and variances. This corre- 
sponds to expanding in the small error from decoherence 
at short times and ignoring the curvature of the Bloch 
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sphere for sufficiently short times, when the spin uncer- 
tainty distribution remains on a locally flat region of the 
Bloch sphere. To linearize the equations we assume that 
all (connected) correlations of order higher than two van- 
ish. The linearized equations are valid for short times, so 
we also make use of the initial conditions in the spin co- 
herent state at t = 0, which are (J y ) = (J z ) = (C yz ) = 



and (jy) = (J 2 ) = J/2. Here we define the covari- 
ance operator C yz = (J y J z + J z J y )/2, while in the main 
text we refer only to its average, V yz = (C yz ). Within 
these approximations, and using Eqs. (l25|) and(27), the 



linearized equations for the spin averages required to cal- 
culate the squeezing parameter £ 2 are 



d t (J x ) = -r 2 (J x ) (28) 

d t (J y ) = xj (j z ) - r 2 (j y ) - r 7 (n th + i) (j y ) + r 7 (c yz ) (29) 

d t (J z ) = -2r T (n th + 1) (J z ) - r 7 [J(J + 1) - (J 2 Z )} (30) 

d t (J 2 y ) = 2JA (C yz ) - 2T 2 J 2 ) - - 2r 7 (n th + §) (J 2 - J 2 ) + T 7 J (J z ) + ^ (J.) (31) 

fit < Jl) = -2r T (n th + [3 ( J 2 Z ) - J(J + 1)] + T 7 (J z ) [l - 2J(J + i)] (32) 

d t (c yz ) = xj (Jl) - r 2 (c yz ) - sr 7 (n tu + ±) (c yz ) - r 7 (j 2 - \) (j v ) (33) 



This linear set of equations can be directly solved. The 
full analytic solutions are lengthy so we simply plot the 
numerical solution for the squeezing parameter (blue 
dash-dotted in Fig. [4]) , which agrees with exact numerics 
at short times. 

We use these linearized equations to estimate the scal- 
ing of the optimal squeezing parameter (see Eqs. (7) and 
(8) in the main text). First, we solve Eqs. p8p3] ) to 
second order in t. Second, we calculate £ 2 (see Eq. (5) in 
the main text) from the resulting spin averages, and sim- 
plify the result in the limit of interest, J 3> l,«th ^ 1, 
and assuming JXt ^> 1 as required for significant squeez- 
ing. Third, we assume that all sources of decoherence are 
small, and expand in the errors T^n^t <C 1 and T 2 t <C 1. 
Within these approximations we obtain 



f + 4jr 7 n th £ 
(JXt) 2 



+ (5r 7 n th +T 2 - 1 )t, (34) 



where the first term l/(JXt) 2 is the result from lin- 
earized unitary squeezing, and the remaining terms are 
the lowest order corrections in both sources of decoher- 
ence. We further approximate T 2 X 3> r 7 n t hj valid for 
sufficiently large detuning A, and JT^n^t 3> f, valid 
self-consistently at the optimal squeezing time and in the 
relevant limit J 3> f). Within these approximations we 
obtain Eq. (7) in the main text. Finally, we optimize £ 2 
with respect to t, obtaining Eq. (8) and t opt given in the 
main text. 



Combining individual dephasing and collective relaxation: 
numerics 

As discussed in the main text, in the absence of a me- 
chanical drive we can neglect phonon number fluctuations 
for a sufficiently large number of tt pulses M. In this case 
the remaining sources of decoherence are intrinsic single 
spin dephasing and collective relaxation induced by me- 
chanical dissipation. These sources of decoherence are 
simple to treat separately but difficult to treat simulta- 
neously for a large number of spins. To calculate the 
solid blue curves in Fig. 2 of the main text, we treat 
the combination of both sources of decoherence as ap- 
proximately independent, valid provided both are small 
enough to still allow spin squeezing. To calculate the spin 
averages needed for the squeezing parameter, we first ac- 
count for collective relaxation using the Dicke state basis 
in which total J is conserved. We then account for in- 
dividual dephasing by multiplying the resulting averages 
by dephasing factors such as (J x (t)) = e~*/ T2 (J x (t))n 
where (Jx(t))r> 1S the result of the Dicke state calcula- 
tion. Each step would be numerically exact in the ab- 
sence of the other source of decoherence; thus we expect 
that this procedure provides a good approximation if all 
errors are small. To verify the accuracy of the approach, 
we compare the result with exact numerics calculated by 
numerically integrating the full master equation for small 
N in Fig. g] 



Phonon number fluctutions 



In this section we consider fluctuations of the phonon 
number, n — a^a. We start by rewriting the effective 




FIG. 4. Spin squeezing parameter versus scaled precession 
time with N = 8 spins. Solid red lines show squeezing calcu- 
lated using the approximation discussed in the text, treating 
single spin dephasing and collective relaxation independently, 
with nth = 0,10,100 (bottom to top). Red dots show the 
exact numerics. The detuning A is optimized for each value 
of nth. Blue dash-dotted lines show squeezing from linearized 
equations. Solid black line shows unitary squeezing, dashed 
black line shows unitary squeezing from linearized equations. 
Parameters are uj m /2n = 1 GHz, gjlix = 1 kHz, Q — 10 6 , 
T 2 = 100 ms. 



Hamiltonian [see Eq. (4) in main text] for the collective 
spin coupled to a driven oscillator in the frame rotating 
at the mechanical drive frequency, 



H cS = AJ Z + \J z a)a + ^J 2 



- 5a)a + £l{a + a}) (35) 



where A = Ab — <jJd is the detuning of the magnetic tran- 
sition frequency from the drive, and 5 = Ud — u> m is the 
drive detuning from the mechanical frequency. The am- 
plitude of the drive is fi and we have made the rotating 
wave approximation. Our aim is to find the effect of the 
oscillator on the spin to second order in A (within the 
Gaussian approximation). For this we require the num- 
ber fluctuation spectrum of a damped, driven, thermal 
oscillator in the absence of coupling to the spin. 



Number fluctuations of a driven thermal mode 

To calculate the effective dephasing from number fluc- 
tuations, we first need the power spectral density of 
phonon number fluctuations, 



S n {u) = / dte luJt (6n{t)Sn{0)) 



(36) 



where 



i^a, Sn = 



(n) and the average is taken 



with respect to the oscillator in thermal equilibrium with 
its environment. In the absence of coupling, A = 0, the 
Langevin equation for the driven thermal mode in the 
frame of the classical drive frequency and within the ro- 
tating wave approximation is 



a(t) = (iS - 1) a(t) + n + 



(37) 



The solution is a(t) — a + d(t), where a = _ iS +^/ 2 
coherent amplitude due to the drive, and 



is the 



(38) 



describes thermal and quantum fluctuations. The mean 
phonon number is the sum of driven and thermal parts, 
n = njr + n-th, where 



O 2 



n dr 



a 



8 2 + 7 2 /4 



(39) 



and the thermal occupation is n t h = (d^d^ = l/(e Wm / T - 
1). Using Eq. ([38| we find the two-time correlations, 



{S{t)d{Q))=f hb e- iSt e-^'\ 
(d{t)d\Q)) = {n th + l)e m e-^l 2 



(40) 
(41) 



and from these we can calculate the full spectrum of 
driven thermal number fluctuations. The correlation us- 
ing Wick's theorem and a(t) = a + d(t) is 

((Sn{t)Sn(Q))) = (n(t)n(O)) - n 2 

=n dr [(rf(t)d(0)) + (d(t)d\0))] (42) 
+ (d\t)d(0)) (d(t)d\0)) . (43) 



Using Eq. ( 38 1 and taking the Fourier transform, we find 



that the number fluctuation spectrum for a driven ther- 
mal oscillator is given by 



nth »th + 1 

_(u;-<5) 2 + 7 2 /4 (lj + 6) 2 

27"th(^th + 1) 



lu 2 + 7 2 



(44) 



Effects of number fluctuations on the spin in the Gaussian 
approximation 

From the spectrum of number fluctuations we can cal- 
culate the effect of number fluctutions on the spin de- 
phasing and squeezing. We write the full Hamiltonian as 
H c g = Hq + H osc + H lnt , where H osc describes the driven 
damped oscillator, and 



(45) 



describes the spin including the constant effective squeez- 
ing term. The coupling in the interaction picture and in 
the toggling frame is 



flint (*) = Xf(t)Sn(t)J z , 



(46) 



where J z is time-independent as it commutes with the 
full Hamiltonian. We have included the function f(t) to 
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describe spin echo, which effectively inverts the sign of 
the interaction with each 7r pulse. 

The equation of motion for the operator J + in the in- 
teraction picture is 



j+(t)=i\f(t)[J z 6n(t),J + (t)]. 



(47) 



We integrate this formally, insert the solution, and take 
the average with respect to the oscillator to get 



J + (t) = -\ 2 f(t) I f(t')([j z Sn(t)[JJn(t'),J + (t')]}) osc 

(48) 



where (-) osc denotes averaging over the oscillator degrees 
of freedom. Note that we neglected additional noise 
terms; these play no role as we will only be interested in 
taking the average at the end. Next, we neglect the time 
dependence of J+{t') under the integral, as it is higher 
order in A, J+(i) = e lHot J + e~ iHot = e i ^ J ^ 1 ^ t J + . Ex- 
panding the commutators we obtain 



= -y.f(t) / f(t')[J z J + (5n(t)Sn(t')) - J + J z (5n(t')Sn(t))] 
Jo 

= -\ 2 f{t) f f(t') f ^e-^-V [J,J + S n (v) - J+JM-u)] 



using Eq. (36). Defining the symmetric and antisymmetric parts of the number fluctuation spectrum, 



S n (uj) = - [S n (uj) + S n (-u))] , A n (oj) = S n (oj) - S n (-uj), 



we can rewrite Eq. (49 1 as 



j+ = -\ 2 f(t)J*f(t>) [ ' L ', < v 



2tt 



S n (u>) + A n (u>) J z - 



1 



Solving Eq. (51 1 and finally taking the average with respect to spin degrees of freedom, we obtain 

(J + (i))=e-*(e^-*)j + (0)) 
where (•) is the average over all degrees of freedom, and 

X-A 2 | J dt x j dtte-toto-^Ofa-tJfitMfo), 



(49) 

(50) 

(51) 

(52) 

(53) 
(54) 



Here all integration limits are from — oo to oo, and 



To evaluate x an d A* f° r a given pulse sequence, we 



the time integration limits are accounted for in fit') oc next define /(w) = J dte lU3t f(t) to rewrite the double 

6(t')6(t—t') and the step function 0(ti — t 2 )- Similarly, we time integral as 
obtain the other averages needed to calculate the squeez- 
ing, /, 



(jl(t)) = e-^(e^-Vjl(0)), (55) 
(J+(t)J z (t)) = e-x (e^( J «"4) J+(0)J,(0)) . (56) 

By comparing with the spin evolution under unitary one- 
axis twisting, we see that fx describes spin squeezing with 
an effective squeezing coefficient A c ff = fj,/t. The param- 
eter x describes collective dephasing. 



dh J dhe-^-^eih-t^fih jfiU) (57) 
duj\ 



2n 



i/K)r 



tt8{w 4- uj\) 



UJ + LOl 



(58) 



We define the filter function for pulse sequence with time 
r between ir pulses, 



F(ct) = ^|/H| 2 . 



(59) 
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The dephasing term is 



Dephasing and squeezing from driven thermal oscillator 



0\2 / duJ a I \ f dw l F ( W 1 T ) 

X = 2A / —S n (uj) / 

2n J 2tt lo( 



TrS^LO + u>i) 



LO + LO\ 

(60) 



Since F{lot) and S n (uj) are both even in u>, the imaginary 
part of the integrand is odd and integrates to zero. As a 
result x is rea l an d we obtain 



X = A 
The coherent term is 



dto F(lot) 

^TT^O 2 



2 < "~ ' V "^5„(w). 



(61) 



fi = i\ / — A n (fjj) 



■k5(lo + LO\) 



to + u>i 
(62) 



Since j4 n (w) is odd, in this case the imaginary part (in- 
volving the (5-function) is zero. The real part is 



, 2 f duo K(lot) . 



where we defined 

K[ut) = 2lj 



2 I du>i F{w\t) 
2tt Luf(oj + LUl) 



(63) 



(64) 



K(ujt) and F(lut) satisfy a Kramers-Kronig relation 
(with a factor of u> 2 from the definitions). 



Adding a mechanical drive, the total dephasing from 
number fluctuations becomes \ — Xth + Xdr, where %dr is 
obtained from the driven part of the number fluctuations, 



(w-5) 2 + 7 2 /4 (w + (5) 2 + 7 2 /4 
The dephasing involves the symmetrized part 

1 



(67) 



S n T {ui) =7TCdr Thb. + 



1 



1 



(w- ( 5) 2 + 7 2 /4 (lo + S) 2 



(68) 



Using Eq. (61) this yields the dephasing from a driven 



for an M pulse sequence. In the relevant limit 71 < 1 
and choosing the timing r = 2tt/5, wc obtain Xdr given 
in the main text. 

For a driven oscillator the power spectral density is 
not symmetric, and the asymmetric part can lead to ad- 
ditional squeezing. The asymmetric part of S n (ui) is 



A n (co) 



7"dr 



1 



(cj + 5) 2 +7 2 /4 (w-<5) 2 +7 2 /4 



(69) 

Using Eq. ( 63 ) , and choosing the pulse timing to coincide 



with a coherence "revival", r = 2tt/S, and assuming the 
mechanical Q > 1, we obtain /id r given in the main text. 



Optimized squeezing with drive 



Dephasing from purely thermal oscillator 



With a strong mechanical drive, the approximate opti- 
mal squeezing is obtained similarly as in Sec. [] above. In 
the driven case we assume that the detuning is large, so 
From Eq. (|44j), the number fluctuation spectrum of a that T 7 -> 0, and the drive is strong so that n dr > n th 

and Xdr 3> Xth- Again expanding £ 2 in the limit J 3> 1, 
^th 3> 1, and small errors tjT-x-, Xdr < 1, we obtain 



purely thermal oscillator in the frame of the mechanical 
drive is 



27n th (n t h + 1) 



■T 



(65) 



The spectrum is symmetric in frequency and thus /i = 0. 
We obtain the dephasing from thermal flucuations by 



inserting Eq. (65) in Eq. (61). For a sequence with an 



even number of pulses M, the filter function is 



7 n th 



t 



(70) 



where we chose Udj ~ (S/g) 2 , the maximum allowed driv- 
ing strength in our perturbative treatment of the cou- 
pling. Optimizing with respect to t we recover Eq. (8) in 
the main text. 



Fm {lot) = 2 sin f — ^— J 



(LOT 

1 — sec — 
V 2 



(66) 



where r is the time between evenly spaced pulses and 
the total sequence time is t = Nr. In the relevant limit 
7T < 1 we obtain Xth given in the main text. 
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